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Packing is a complex phenomenon of prominence in many natural and industrial processes (liquid crystals, 
granular materials, infiltration, melting, flow, sintering, segregation, sedimentation, compaction, etc.). A variety 
of computational methods is available in particular for spheroid particles. Our aim is to apply the principle 
of the random sequential addition algorithm but with small blocks of varying size and orientation. Here the 
main purpose is to reproduce the observed arrangement of graphitic assemblies in polymeric matrices. Random 
packing is improved by applying an external pressure implemented with a drifted diffusive motion of the fillers. 
Attention is also paid to the emergence of structural and orientational order. Interestingly, mixtures of fillers of 
irregular shapes can be dealt with efficiently using the proposed algorithm. 
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I. INTRODUCTION 

In the present Article we are faced with the problem of 
realizing the densest packing with small blocks out of a bi- 
modal size distribution. The motivation is to generate the ar- 
rangement of graphitic assemblies in polymeric matrices as 
obtained after compaction (see Fig. [TJ. In these materials, the 
graphitic phase can reach high volume fractions. It follows 
that the packing structures of these systems demand attention, 
especially when the electrical properties are under investiga- 
tion, since the polymer is insulating [1,2]. More generally, the 
polymeric matrix is enriched with a second dispersed phase in 
order to enhance a wide range of properties of the resulting 
composite material (mechanical, electrical and thermal) Q. 
In the appearance, this kind of materials display quite disor- 
dered microstructures that nonetheless reveal attributes of key 
importance in several aspects J4). At the extremes, random 
packing can be addressed with detailed, fully deterministic 
numerical methods like classical molecular dynamics [5 6| 
or with fully random processes QUI. In between, dynamic 
approaches allowing the particles to change their size 191 4121 
appear to be the most established. Typically, it arises that the 
shape of the fillers has a significant influence on the packing 
properties (see in particular Refs. J9][T3]). Unfortunately, the 
mathematical treatment of general morphologies is far from 
being trivial and for this reason most of the research has so 
far concentrated on ellipsoidal particles (spheres, oblates and 
prolates), and in a few exceptions on cubes lfT4l4T6ll and oth- 
ers ifTTll . Our work is based on the random sequential addition 
(RSA) algorithm |7| with the possibility to use hard fillers of 
arbitrary shape, similarly to what reported in Refs. lfT8Tl20l 
for 2D, but with coarser discretization for efficiency concerns. 
In order to favor compaction, we introduce the effect of an 
external pressure via a biased diffusive motion of the fillers. 
In principle, blocks can pack very well because they have ex- 
tended planar surfaces [9|. When they can be rotated, severe 
geometrical constraints are imposed in their neighborhood. 
For this reason, the structural properties of the final config- 
urations are analyzed in detail by means of radial distribution 
functions. Furthermore, we work out a method to detect ori- 
entational correlations at a local level. Our results suggest that 



orientational correlations have a range that can go beyond that 
of structural order. 



II. RANDOM ADDITION STRATEGY 

The domain to be filled with small, impenetrable blocks has 
all side lengths equal to 1 with periodic boundary conditions. 
We build the blocks from small spheres of radius r, for the 
reasons explained in Ref. 1 18 1. The blocks are obtained as fol- 
lows. A single sphere is first placed randomly in the domain 
with uniform distribution. This sphere is that with the smallest 
x, y and z coordinates composing the block (before rotation). 
The other spheres are added so as to be aligned along the three 
cartesian axes with their centers separated by the distance r. 
Let l x , l y and L be the side lengths of the block; for example, 
the length of the x side is l x = (n x + l)r, where n x is the number 
of spheres aligned along this direction. The size of the block 
is fixed randomly by choosing n x , n y and n z from three Gauss 




Figure 1 : Microstructure geometry of a typical polymeric matrix re- 
inforced with a graphitic powder. It appears that the graphitic phase 
(lighter) displays to a certain extent rectangular characteristics. By 
2D image processing it is determined a volume fraction of around 
62%. 
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Figure 2: How the spheres are inscribed in the cubes. The origin of 
axes is in the center of the first placed sphere of the block. In red is 
shown the boundary of the block in the xy plane in the absence of tilt. 



distributions of means (n x ) = ((l x ) —r)/r, (n y ) and (n z ); the 
three distributions have the same variance a 2 and subsequent 
computations are done using the integer parts of n x , n y and n z . 
The resulting block is then rotated around the z, x and y axes, 
in this order, by the angles 9 Z , 9 X and 9 y , respectively. These 
three angles of rotation are chosen randomly within a given 
interval from a uniform distribution. For the rotation, the axes 
are translated with their origin in the center of the first sphere. 

Let Ax, Ay and Az be the distances in absolute value along 
the three axes between the centers of two spheres. This means 
that the two centers are separated by yj (Ax) 2 + (Ay) 2 + (Az) 2 . 
The block is added to the domain if for every center of its 
spheres at least one A is larger than 2r + d m [ n [21 1. If this con- 
dition is not fulfilled there is overlap: the block is discarded 
and the above procedure is repeated. 

We want to verify relatively fast the overlap condition. To 
this end, the domain is divided into five layers of equal width 
parallel to the xy plane. The centers of the spheres are stored 
in five lists, one for every layer. In order to avoid overlap at 
the interfaces, the spheres falling within a distance of 2r from 
an interface are also added to the list of the adjacent layer. In 
this way, the overlap condition is checked only among spheres 
belonging to the same list. 

The binning of the domain is regular, leading to A^ 3 voxels. 
Every side of the domain is thus divided into N bins. Know- 
ing the centers of the spheres, the real blocks are readily ob- 
tained at the end of the simulation by inscribing every sphere 
in a cube of side length 2r. Every voxel inside the cubes will 
be considered as occupied by the filler phase. Voxels out- 
side the cube will be considered instead as occupied by the 
matrix phase. By this procedure, the surface of the resulting 
parallelepipeds is rough. The degree of roughness is arguably 
in line with experimental observations. Indeed, we find that 
the depth of the defects is given approximately by the for- 



mula r | cos 9 Z sin 9 Z |. The roughness is of course maximal for 
9 Z = k/A. By simple geometrical arguments, we find that the 
depth of the incidents at the surface is approximately at most 
r/2 (see Fig. |2j». This value corresponds to 1% of the side of 
the simulation domain for r = 0.02. Indicatively, the side of 
the simulation domain can be assumed to be 350 [firm], lead- 
ing to a depth of the defects of the surface of 3.5 [fim]. In the 
case of graphite-polymer composites, we observe that at the 
micron scale the surface of the fillers is irregular (see Fig. [TJ. 

Diffusive motion is introduced by means of random moves 
of the blocks along the x, y and z directions. The effect of an 
external pressure is accounted for by introducing a bias toward 
the center of the domain. More precisely, after a given num- 
ber of unsuccessful attempts, all the already placed blocks are 
moved randomly. For every block, first the direction of the 
trial move is chosen randomly. The sense of the move is then 
chosen with probabilities p x , p y and p z . For example, for a 
block with the first sphere in the region defined by x > 1/2, 
p x is the probability to be translated toward the center with 
sense and direction defined by the vector (—1,0,0). Given a 
block, the distance spanned by the trial moves is different for 
every direction. When a move is not possible because result- 
ing in overlap, the distance associated with that direction is 
reduced by a factor of 2. The first time one of the possible 
move distances falls below <i m i n , the block is no more taken 
into account for random moves. Finally, in this way, purely 
diffusive motion corresponds to p x = p y — p z = 0.5 while in 
all other cases there appears a resultant pressure (inward or 
outward oriented). 



III. SIMULATION SETTINGS AND MODELING 
APPROACH 

In all simulations reported here, the standard deviation a is 
set to the value y/03 and the radius of the spheres is r — 0.02. 
We always start with a given average size and smaller blocks 
are considered when no block can be placed after 200'000 
trials. The simulation is terminated when no smaller block 
can be placed for the same condition. Table [I] lists the pow- 
ders with the average size of blocks they are composed of. 
Mixtures and aspect ratios are chosen so as to obtain com- 
posite materials produced experimentally by using common 
graphitic powders. 

Different random orientations are investigated. We distin- 
guish the types of arrangement listed in Tab. [II] Every system 
is considered without diffusion and pressure and with pres- 
sure. In the latter case, for the probabilities of rearrangement 
we choose p x — p y = p z = 1 in order to reach the maximal 
packing fractions. The initial possible displacement along the 
three directions is 0.5. The minimal distance of approach is 
set to c/min = 2/N, with N = 256, corresponding to the width 
of two adjacent voxels. The blocks are never moved twice 
consecutively, but one unsuccessful attempt is sufficient in or- 
der to call the routine performing diffusive motion. The aim 
is to keep the pressure effect as close as possible to the lin- 
ear dependence on the total number of blocks when present in 
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Particle 1 Particle2 


Volume ratio 


bimodall 


(0.2,0.2,0.12) (0.06,0.06,0.04) 


33.3 : 1 


bimodal2 


(0.2,0.2,0.12) (0.1,0.1,0.06) 


8 : 1 


bimodaB 


(0.2,0.2,0.12) (0.14,0.14,0.08) 


3.1 : 1 


granular 


(0.2,0.2,0.12) (0.2,0.2,0.12) 


1 : 1 


fine 


(0.06, 0.06, 0.04) (0.06, 0.06, 0.04) 


1 : 1 



Table I: Generated powders of different mixtures. The vectors (/ t , / v , / z ) defines the average side lengths of the blocks. 



No. 


Type 


0, 




V 


6; 


I 


aligned 













II 


misaligned 










[-rc/4,+;r/4] 


III 


tilted x 


[-7T/36H 


-7t/36] 








IV 


tilted xy 


[-7T/36H 


-7t/36] 


-ff/36,+7E/36] 





V 


misaligned and tilted x 


[-7T/36H 


-7t/36] 





[-rc/4,+;r/4] 


VI 


misaligned and tilted xy 


[-n/36H 


-tc/36] 


-tc/36,+w/36] 


[-^/4,+^/4] 



Table II: Case studies and their orientation characteristics. 





I II III IV V VI 


bimodall (std) 


55.1% 49.5% 54.4% 52.5% 50.1% 49.4% 


bimodall (press) 


59.5% 54.3% 57.4% 57.2% 54.9% 52.1% 


bimodal2 (std) 


52.0% 49.3% 52.2% 51.0% 47.2% 47.9% 


bimodal2 (press) 


58.7% 51.4% 56.7% 55.7% 51.0% 50.3% 


bimodaB (std) 


50.4% 45.0% 48.8% 48.8% 43.1% 42.9% 


bimodaB (press) 


54.9% 47.6% 52.9% 51.5% 46.1% 46.5% 


granular (std) 


40.9% 37.0% 38.8% 41.5% 35.7% 34.7% 


granular (press) 


44.5% 40.3% 43.0% 41.3% 42.4% 37.9% 


fine (std) 


34.1% 33.4% 33.8% 33.8% 32.8% 33.3% 


fine (press) 


36.6% 35.7% 36.4% 36.3% 35.8% 35.7% 



Table III: Packing fractions as obtained from simulations for systems with different orientational characteristics and size distributions. The 
name of the mixtures is defined in Tab.|l] The labels I- VI refer to the systems with the degree of orientational randomness defined in Tab, [ill 
The abbreviation std stands for standard (without pressure and diffusive motion); press means that an external pressure is applied as explained 
inSec.|IIl] 




Figure 3: Sections of some packing structures for the mixture bimodall (see Tab. [TJ with an applied external pressure, a) Aligned case, b) 
Misaligned case, c) Tilted xy case. Note that the boundary of the blocks in the plane xy can be diagonal, d) Misaligned and tilted xy case. Note 
that the boundary of the blocks can be irregular in this case. 
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I 

no. part, density 


II 

no. part, density 


III 

no. part, density 


IV 

no. part, density 


V 

no. part, density 


VI 

no. part, density 


, , is large 
bimodall (std) 5 

small 


101 41.8% 
1'222 13.3% 


81 34.2% 
1'423 15.3% 


93 40.1% 
1'325 14.3% 


91 37.9% 
1'345 14.6% 


81 35.0% 
1'436 15.1% 


77 34.2% 
1'448 15.2% 


bimodall (press) ^ ar 8 e 
small 


110 45.4% 
1'263 14.1% 


88 38.9% 
1'405 15.4% 


100 41.8% 
1'415 15.6% 


102 42.9% 
1'290 14.3% 


92 41.1% 
1'271 13.8% 


81 35.9% 
1'498 16.2% 


bimodal2 (std) ge 
small 


93 38.1% 
346 13.9% 


85 36.7% 
310 12.6% 


90 37.7% 
368 14.5% 


90 37.4% 
335 13.6% 


81 34.6% 
288 12.6% 


83 36.0% 
292 11.9% 


bimodal2 (press) ^ ar 8 e 
small 


110 45.7% 
325 13.0% 


88 38.5% 
310 12.9% 


103 43.4% 
334 13.3% 


100 43.2% 
299 12.5% 


91 40.0% 
265 11.0% 


88 38.9% 
275 11.4% 


bimodaB (std) ge 
small 


98 41.1% 
80 9.3% 


85 36.9% 
65 8.1% 


88 37.7% 
96 11.1% 


89 37.9% 
93 10.9% 


81 35.9% 
57 7.2% 


76 33.7% 
73 9.2% 


bimodaB (press) ^ ar S e 
small 


104 43.5% 
98 11.4% 


89 38.6% 
71 9.0% 


98 42.1% 
92 10.8% 


103 44.1% 
61 7.4% 


92 40.2% 
47 5.9% 


84 38.3% 
66 8.2% 


, ... large 
granular (std) 

small 


99 40.9% 
0% 


88 37.0% 
0% 


92 38.8% 
0% 


98 41.5% 

0% 


81 35.7% 
0% 


80 34.7% 
0% 


granular (press) ^ ar 8 e 
small 


107 44.5% 
0% 


93 40.3% 
0% 


104 43.0% 
0% 


98 41.3% 
0% 


90 42.4% 
0% 


86 37.9% 
0% 


fine (std) 

small 


0% 
3'066 34.1% 


0% 
3'000 33.4% 


0% 
3'014 33.8% 


0% 
3'063 33.8% 


0% 
2'962 32.8% 


0% 
2'993 33.3% 


fine (press) large 
small 


0% 
3'231 36.6% 


0% 
3' 148 35.7% 


0% 
3' 182 36.4% 


0% 
3' 191 36.3% 


0% 
3'144 35.8% 




3'138 35.7% 



Table IV: Number of block particles of type larger and smaller composing every bimodal powder and their contribution to the packing density. 
Misalignment of larger blocks is the primary cause of significant variations of the overall volume fraction. The mixtures and the size of their 
particles are introduced in Tab.|l] Roman numbers specify the degree of orientational randomness as defined in Tab.|TI] 



relatively large number. 

Of course, the application of an external load has the ef- 
fect to deform the fillers. Furthermore, anisotropic pressures 
can also lead to a certain degree of orientational order. We 
apply an isotropic pressure, but the shape of the fillers and 
the restricted range of the tilt angles intend to simulate the ef- 
fect of uniaxial compaction along the z direction. Graphitic 
powders are the result of an extremely expensive process of 
milling. Because of their molecular structure, the size of par- 
ticles is reduced through both exfoliation and fragmentation. 
It is clear that the combination of these phenomena introduces 
faces and angularities in the resulting particles. In that respect, 
the block is the simplest geometry approaching the shape of 
real graphitic assemblies. 



IV. RESULTS AND DISCUSSION 




In Tab. Ill the packing fractions <j) reached for the different 
powders and orientational characteristics are reported. In gen- 
eral, the application of pressure leads to denser packings. We 
ascribe the unique exception to statistical fluctuations. Indeed, 
for the granular powder of type IV with pressure, average over 
five realizations gives <p = 0.422, in this case above the value 
for the same system with no applied pressure. For compara- 
tive purposes relative to the different orientational properties, 



Figure 4: 3D rendering of the final structure for the powder bimodal2 
with pressure in the aligned case. The volume fraction is 58.7%. 



Fig. [3] shows the sections of four final configurations for the 
powder bimodall with an applied pressure. In this case, the 
highest packing values are reached since the smaller particles 
can better fill the irregular interstices between larger blocks. 
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I 


II 


III 


IV 


V 


VI 




0-0 


0.5159-0.1329 


0-0 


0-0 


0.5209-0.1355 


0.5072-0.1287 


bimodall (std) 


0-0 


0-0 


0.0586-0.0017 


0.0588-0.0017 


0.0585-0.0017 


0.0590-0.0017 




0-0 


0-0 


0-0 


0.0576-0.0017 


0-0 


0.0585-0.0017 




0-0 


0.5140-0.1321 


0-0 


0-0 


0.5142-0.1324 


0.5084-0.1294 


bimodall (press) 


0-0 


0-0 


0.0579-0.0017 


0.0577-0.0017 


0.0570-0.0016 


0.0585-0.0017 




0-0 


0-0 


0-0 


0.0581-0.0017 


0-0 


0.0580-0.0017 




0-0 


0.4895-0.1206 


0-0 


0-0 


0.4761-0.1153 


0.4808-0.1182 


bimodal2 (std) 


0-0 


0-0 


0.0571-0.0016 


0.0565-0.0016 


0.0550-0.0015 


0.0564-0.0016 




0-0 


0-0 


0-0 


0.0553-0.0015 


0-0 


0.0566-0.0016 




0-0 


0.4780-0.1146 


0-0 


0-0 


0.4884-0.1198 


0.4935-0.1218 


bimodal2 (press) 


0-0 


0-0 


0.0584-0.0017 


0.0586-0.0017 


0.0553-0.0015 


0.0553-0.0015 




0-0 


0-0 


0-0 


0.0585-0.0017 


0-0 


0.0565-0.0016 




0-0 


0.5001-0.1248 


0-0 


0-0 


0.4896-0.1227 


0.4707-0.1114 


bimodaB (std) 


0-0 


0-0 


0.0579-0.0017 


0.0533-0.0014 


0.0493-0.0012 


0.0568-0.0016 




0-0 


0-0 


0-0 


0.0570-0.0016 


0-0 


0.0583-0.0017 




0-0 


0.4870-0.1186 


0-0 


0-0 


0.4694-0.1127 


0.4577-0.1077 


bimodaB (press) 


0-0 


0-0 


0.0545-0.0015 


0.0566-0.0016 


0.0538-0.0016 


0.0590-0.0017 




0-0 


0-0 


0-0 


0.0524-0.0014 


0-0 


0.0556-0.0015 




0-0 


0.4949-0.1200 


0-0 


0-0 


0.4870-0.1222 


0.4779-0.1185 


granular (std) 


0-0 


0-0 


0.0555-0.0015 


0.0571-0.0016 


0.0531-0.0014 


0.0532-0.0014 




0-0 


0-0 


0-0 


0.0548-0.0015 


0-0 


0.0579-0.0016 




0-0 


0.5175-0.1336 


0-0 


0-0 


0.4617-0.1126 


0.5162-0.1312 


granular (press) 


0-0 


0-0 


0.0579-0.0016 


0.0537-0.0014 


0.0609-0.0018 


0.0601-0.0018 




0-0 


0-0 


0-0 


0.0577-0.0016 


0-0 


0.0570-0.0016 




0-0 


0.5227-0.1366 


0-0 


0-0 


0.5162-0.1332 


0.5146-0.1325 


fine (std) 


0-0 


0-0 


0.0589-0.0017 


0.0581-0.0017 


0.0581-0.0017 


0.0581-0.0017 




0-0 


0-0 


0-0 


0.0582-0.0017 


0-0 


0.0576-0.0017 




0-0 


0.5200-0.1351 


0-0 


0-0 


0.5162-0.1331 


0.5095-0.1299 


fine (press) 


0-0 


0-0 


0.0589-0.0017 


0.0581-0.0017 


0.0578-0.0017 


0.0587-0.0017 




0-0 


0-0 


0-0 


0.0574-0.0016 


0-0 


0.0587-0.0017 



Table V: Average value and variance of the variables AO for the case studies of Tabs. |l| and [TT] For every mixture, the first row refers to AO- 
followed by its variance, the second to A0 X and its variance, and the third to A0 y and its variance. The abbreviations std and press mean with 
no pressure and with pressure, respectively. 



Figure [4] shows the 3D representation of the packing structure 
when the role of pressure is more effective (bimodal2 powder 
of type I). In general, orientational disorder lowers the packing 
density in particular if misalignment is allowed (see Fig. BJ. 
Fine powders lead to looser random packings. It turns out 
that the highest volume fractions can be achieved with more 
polydisperse size distributions by starting from the addition 
of larger blocks in the absence of orientational disorder. The 
role of smaller particles is to fill progressively narrower inter- 
stices. In this way, it is possible to simulate the infiltration of 
finer powders between stuck particles and eventually their ag- 



gregation. Table IV further details the results for the packing 



structures of all powders. 

It is interesting to look at how the packing fraction pro- 
gresses in the course of time. Figure [6] shows the cases of 



granular and fine powders of type I (aligned case). With pres- 
sure, the volume fraction follows a power law (j) ~ t V over one 
more decade. Furthermore, at the end, these simulations are 
faster. For other orientational characteristics, it is found that (j) 
displays an analogous behavior. In the case of bimodal pow- 
ders, we also find similar curves for the addition of smaller 
blocks. 

In order to gain insight into the structural properties of the 
final configurations of the generated powders, we calculate the 
radial distribution functions (RDFs) between types of blocks. 
The usual definition of the RDF g(r) is [22 1 



-J-g(r)4nr 2 Ar = S(r) 



(1) 



V = 1 is the volume of the domain. S(r) is the average number 
of particles of type 2 falling within a spherical shell of radius 
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Figure 5: Packing fraction against the volume ratio; pressure is ap- 
plied. Aligned case (x), misaligned case (+) and misaligned with 
tilt xy case (*) (cf. Fig.[3jl. Inset: A(j) = 0p re ss — 0std as a function of 
the volume ratio. 0p reS s is the packing fraction reached with pressure 
and st( j that realized with no pressure. 




— granular press 

— tine press 

- - -granular sld 

- - - line std 



2 3 



Figure 6: Packing fraction <j) as a function of time (log-log scale). 
All powders are of type aligned (cf. Tab,|n} with pressure (press) and 
with no pressure (std). 



r and width Ar centered around the particles of type 1 present 
in number N\ . (Precisely, the distance between two blocks is 
determined using their centers of mass.) As a consequence, 
the RDF is a measure of the average radial dependence of the 
density of particles of type 2 around those of type 1 . In the fol- 
lowing, given a powder, with the notation gi s (r) we intend the 
RDF of smaller particles around larger ones. We employ sim- 
ilar notations for the other cases and, with no upper-script, the 
RDF is among all particles. Figure [7] shows different RDFs 
in several cases. In the first row, the main peak of g\\ is de- 
termined by adjacent in-plane blocks. The first peak around 
r = 0.1 is due to the elongated shape of the blocks. It is lower 
than the second peak because a block can have indicatively 
twice more in-plane neighbors. In the absence of pressure 
(second row), it is seen that the peaks are lower and slightly 
broader, indicating that these structures are less compact. Of 
course, these remarks hold also for the other bimodal pow- 



ders. The third row of Fig.|7]shows g\ s for different mixtures. 
Mainly, as the volume ratio gets smaller, there appears more 
distinctly a peak at r « 0. 1 . This means that the contribution of 
out-of-plane small particles become statistically more signifi- 
cant. For the powder bimodaO, the last peak around r = 0.4 is 
related to the first peak of g\\ and corresponds to pairs of large- 
small blocks separated by a large one. For the bimodall pow- 
der, the first peak of g ss always occurs before r — 0. 1 (fourth 
row). It follows that this peak is due to the small particles fill- 
ing the voids between larger blocks. The same RDF for the 
fine powder is similar (fifth row). The absence of interstices 
between larger blocks has no consequence because the small 
particles tend to fill densely the voids. Indeed, for the other 
bimodal powders there appear additional peaks after r = 0.1, 
more distinctly for the lowest volume ratio (bimodal3 pow- 
der). The RDF among all particles for the bimodall powder 
(last row) tells us that the statistics is dominated by the small 
particles relatively close to each other determining the first 
peak. In the other cases, the arrangement of smaller blocks 
is progressively similar to that of larger ones and there appear 
two peaks for r > 0. 1 for the reasons explained above (see first 
row). 

Interestingly, Fig. [3] suggests that the orientation of larger 
blocks might be correlated. In order to better understand this 
point, we carry out a statistics of the rotation angles. Let P 
be the number of possible pairs of blocks. We introduce the 
average quantities 



(A0,)--Lia(/)-0z(;)l; 

i>j 

((A0 z ) 2 ) = ^£[0z(O-0zO-)] 2 - 

i>j 

9 z (i) is the rotation angle of the 2-th block. Similar aver- 
age values are defined for the tilt angles 9 X and 9 y . Besides 
the average values (A9), of interest are also the variances 
((Af9) 2 ) — (A9) 2 . In Tab. |v| these quantities are given in all 
cases addressed by simulations. The blocks are expected to 
have the same orientation to a higher degree when the quanti- 
ties (A9) are small. It is important to remark that (A9 Z ) can be 
at most n/2. It is found that this quantity varies in the interval 
0.4617 - 0.5227, i.e. 26° - 30°. As a result, from this overall 
statistics it clearly arises that orientational ordering is present 
in the internal packing structures. However, it seems difficult 
to discriminate between and compare the different cases on 
the basis of this method. 

The emergence of orientational order can be addressed in a 
more specific way. Let B be the number of blocks. We define 
Bi(r) as the ensemble of blocks inside a spherical shell of ra- 
dius r and width Ar centered around the z'-th block. In order 
to take into account the local, relative orientation of pairs of 
blocks, we introduce the function 

s(r)=y y (r/2)2 

V 7 iT(o^(0-^O , )] 2 + (r/2) 2 ■ 

9 z (i) is the angle of rotation around the z axis of the i-th 
block; the constant F is assumed to be 10~ 4 . In the above 
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Figure 7: Representative RDFs. For Ar we choose 5/256 (see Eq.[TJ. 



sum, the idea is that the contribution of every pair of block 
is "weighted" by a Lorentz function. Summands close to 1 
indicate that the corresponding Z angles of the blocks differ 
slightly. Substitution into Eq. [I] yields a kind of RDF g z (r) 
measuring the average radial dependence of orientational cor- 
relations. Similar definitions hold for the functions g x and g y . 
For completeness, we also consider the function g u defined 
by using the angle formed by the orientation vectors of pairs 
of blocks. For every block, the orientation vector is initially 



u = (1,0,0), to which three rotation transformations of angles 
9 Z , 6 X and 6 y are applied in this order. In Fig.JS] we plot these 
functions for two powders in different cases of orientational 
randomness. For the functions g z of the bimodall powder 
(first row), the peaks around r — 0.3 indicate that the orien- 
tation of every block is to a certain degree correlated with that 
of its adjacent blocks. From the further peaks we can conclude 
that the orientation is also correlated between second-nearest 
neighbors (see Fig. 0J. For the functions g x and g y (second 
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bimodall bimodal2 bimodaO granular fine 



no pressure 


57.9% 


53.2% 


48.3% 


38.5% 


37.9% 


with pressure 


59.8% 


54.8% 


50.3% 


41.6% 


42.0% 



Table VI: Packing fractions reached for powders composed of spher- 
ical particles. The average diameters of the two types of spheres are 
equal to the longest side of the blocks for the corresponding mix- 
ture of Tab. H] Since simulations are much faster in this case, smaller 
spheres are considered after 2 x 10 6 unsuccessful attempts; all other 
simulation settings are maintained. 



row), the average values are always at least a factor 20 larger 
than for the previous g z functions. The reason is that the tilt 
angles are restricted to values belonging to small intervals. 
The presence of more peaks is consistent with the elongated 
shape of particles. For g x all peaks are comparable in height; 
instead, for g y there appears a marked peak around r = 0.3. 
Without pressure (third row), the peaks are reduced in height 
and that near r = 0.3 is always the more pronounced. This 
means that orientational correlations among second-nearest 
neighbors are weaker because of looser packing. For the bi- 
modal3 powder (fourth and fifth rows), from the peaks in the 
range r = 0.1 — 0.2 it turns out that out-of-plane blocks are 
also correlated. This peak is of course higher for the functions 
g x and g y . Furthermore, in this case the functions g x and g y 
under the same conditions are similar (see second row). In 
the last row, the function g u is considered for three different 
powders. Direct comparison with the plots of the functions g x , 
g y and g z for the same type of powders indicate that the first 
peaks flatten out. It follows that, for adjacent in-plane blocks, 
the angles 9 X are weakly correlated; for out-of-plane blocks, 
it arises that the angles 9 Z are instead weakly correlated. 



V. CONCLUSIONS 

It is well-known that the shape and size of the fillers af- 
fect the packing outcome and in turn, for example, the accu- 
rate prediction of the properties of materials by finite element 
methods l27ll . In this work we focused on small blocks be- 
cause they better approximate the shape of graphitic clusters 
in some aspects (edges and corners). The algorithm is based 
on random sequential addition and can be extended easily to 
fillers with more irregular shapes, domains with solid walls 
or random addition mimicking deposition. In this approach, 
the interaction and evolution of the relevant constituents of 
the material are poorly accounted for. The main advantage 
of stochastic methods is the possibility to treat efficiently the 
complexity of the systems by concentrating on simple, ef- 
fective rules (28]. In the case of the present algorithm, full 
control on interparticle distance seems to be the most useful 
property when the potential energy between hard particles is 
assumed to be identically zero, except for the case of con- 
tact. Then, coarser discretization clearly expedites the simu- 
lations. Unavoidably, this leads to defects at the surfaces, but 
roughness is a ubiquitous element in materials science. In the 
case of bimodal packing our results show that more polydis- 
perse size distributions could lead to higher volume fractions, 
because particles are more suited to fill voids comparable to 
their size. In any case, bimodal size distributions deserve at- 
tention for the transport properties associated with their partic- 
ular packing structures. Radial distribution functions provides 
useful information on the size, amount and arrangement of the 
particles composing the systems. Our analysis on correlations 
reveals that structural and orientational order do not coincide. 
It turns out that a nonstructured medium (beyond the radial 
distance r = 0.3 in our systems) can still display forms of or- 
ganization. 



VI. NOMENCLATURE 



We end this section by pointing out that every simulation 
used around 20 MB of RAM memory and the longest run 
took 74h 43min (bimodal3 powder aligned case with pres- 
sure) while the shortest one lasted for 3h 51min (fine pow- 
der tilted x case without pressure). The program is written in 
Perl and it is executed on a desktop computer with an Intel 



processor of the family i7. Finally, Tab. VI provides the pack- 



ing fractions achieved with equivalent spherical fillers. These 
results are comparable with the predictions for random pack- 
ing of spheres ifTTl |23Tl26ll , indicating that volume fractions 
can reach values in the interval 0.60-0.68. It arises that the 
role of pressure is more effective with coarse blocks. The best 
volume fractions can not be reached within the present simu- 
lation settings because we can not tune arbitrarily parameters 
like the aspect ratios, the relative size between small and large 
blocks and between the particles and the simulation domain, 
given the specific application under investigation. It clearly 
turns out that higher volume fractions can be realized only 
with more polydispersity. 



B 

Bi 

^min 
glials >, 
gx,gy,£ 

N 
Ni 

n x ,n y ,r, 



A. Symbols 

number of block particles in the simulation 
domain 

number of blocks in a spherical shell of radius 

r centered around the z'-th block 

minimal distance of approach between two 

block particles 

radial distribution functions 

functions defined for detecting orientational 

order 

side length of block particles 
number of voxels along each side of the sim- 
ulation domain 

number of block particles of type 1 

number of spheres composing each side of 
block particles 
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Figure 8: Functions g x , g y , g z and g u for the mixtures bimodall and bimodaB (cf. Tab.ITl for different orientational characteristics (see Tab.ITTl. 
Ar is set to 5/256 (cf. Eq. [TJ. For the powder bimodall, small particles are omitted. The straight line in red shows the average value. Itis 
calculated by averaging the functions over the bins. 



p number of pairs of blocks in the simulation 

domain 

p x ,Py,Pz - probability for random moves along each axis 

r - radius of spheres composing block particles 

Ar - width of a spherical shell 

g average number of block particles in a spher- 

ical shell of radius r 



t 

u 

x,y,z 
V 



weighted number of block particles in a spher- 
ical shell of radius r 
time 

unit vector 
cartesian coordinates 
volume of simulation domain 
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r - model parameter 

Y - power-law exponent 

- volume fraction 

G - standard deviation of a Gauss distribution 

x ,0 y ,0 z - rotation angles around each axis 



press with pressure 

RDF radial distribution function 

RSA random sequential addition 

std standard 



part. 



B. Abbreviations 



minimum 
particles 
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